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We propose a numerical method for evaluating eigenval- 
ues and eigenfunctions of Schrodinger operators with general 
confining potentials. The method is selective in the sense that 
only the eigenvalue closest to a chosen input energy is found 
through an absolutely-stable relaxation algorithm which has 
rate of convergence infinite. In the case of bistable potentials 
the method allows one to evaluate the fundamental energy 
splitting for a wide range of tunneling rates. 

02.70.Bf, 03.65.Ge 

According to the von Neumann theory [1], in an ideal 
measurement of energy the state of a quantum sys- 
tem collapses instantaneously and completely into some 
eigenstate of the Hamiltonian. If one knows how to 
handle the collapse mechanism and how to select the 
final eigenstate, measurements of energy may be used 
for determining the spectrum of quantum systems. The 
restricted Feynman path-integral approach to quantum 
measurements [2] offers this possibility. During a con- 
tinuous measurement of energy with known result E (we 
consider the particular case of E constant) the Feynman 
paths far from those compatible with the measurement 
result are damped proportionally to the accuracy of the 
measurement itself [3]. Choosing a damping of Gaussian 
type we obtain a Schrodinger problem with an effective 
Hamiltonian 

H eff = H - ihn(H - Ef (1) 

where the k gives the strength of the coupling to the mea- 
surement apparatus. In collaboration with R. Onofrio 
[4] we have recently discussed the dynamics of the wave- 
function collapse induced by the effective Hamiltonian 
(1). Let us consider the case of H with a nondegenerate 
discrete spectrum 

H4> n {x) = E n 4> n {x) (2) 

where in units Ti = 2m = 1 

H = -Vl + V(x) (3) 

with x G R 3 ■ The wavefunction of the measured sys- 
tem can be decomposed in terms of the eigenfunctions 



ip n which are also eigenfunctions of H e f / . Due to the 
presence of the anti-Hermitian term in (1), during the 
measurement the initial wavefunction converges, up to 
a normalization factor, to the eigenfunction with energy 
closest to E at a rate exponentially proportional to k. 

In principle, the numerical solution of the time- 
dependent Schrodinger equation with the effective Hamil- 
tonian (1) represents a relaxation method for evaluating 
eigenvalues and eigenfunctions of H close to the select- 
ing energy E. In practice, we can speed up the relaxation 
by letting k — > oo. In this case we obtain the evolution 
equation 

±iP(x,t;E) = -(H -E) 2 iP(x,t;E) (4) 

where the scaled time t has dimensions [-E 1-2 ]. The func- 
tion ip(x, t; E) is complex in general and we have empha- 
sized its dependence on the selecting energy E. If we 
define the relaxed wavefunction and energy 

iP rel (x;E)= lun rffifjll (5) 
t^oo \\iP{x,t;E)\\ 

E rel {E) = J 4> rel {x-E)*H4> rel {x-E) dx (6) 

they have the property that ip re i(x;E) = ip n (x) and 
E rel {E) = E n when E G T n = }(E n + £„_i)/2, (E n + 
E n+ i)/2[ for n ^ 0. Relaxation to the ground state 
n = is obtained through the weaker condition E G 
r = ]-oc,(£ + £i)/2[. 

It is worth noting the relevance of selectivity. We can 
evaluate whatever eigenstate of the spectrum just giving 
an estimate of it up to an error of the order of the local 
energy spacing. On the other hand, nonselective relax- 
ation methods, like those based on the heat equation [5] 

converge only at the ground state (we suppose H > 0). 
Excited states can be obtained by using an initial wave- 
function orthogonal to all the lower-energy states. How- 
ever, only exactly orthogonal wavefunctions ensure relax- 
ation to the desired eigenfunction. The errors introduced 
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by finite-accuracy numerical orthogonalizations make the 
method unstable and not practical for determining high 
energy states. 

Beside selectivity, another advantage characterizes a 
relaxation method based on Eq. (4). We can solve that 
equation through a finite-difference algorithm which is 
absolutely stable and allows us to evaluate the relaxed 
quantities (5-6) in one step. Let us explain the idea in 
the one dimensional G M. The domain of x can 

be restricted to some interval [x m i„, x max ], depending on 
the selection energy E and the confining potential V(x), 
outside which the relaxed wavefunction vanishes within 
the computer accuracy. The interval is discretized by 
introducing the space lattice 

x Xj = x min + jAx j = 0, 1, 2, . . ., J + 1 (8) 

If the time also is 



where (J + l)Ax = x 
discretized according to 



max ^min • 



t — >■ t m = mAt m = 0,1,2, 



(9) 



Eq. (4) can be reduced to the following set of finite- 
difference equations 



^ ^ = -(7^ + V+^« 1 + 



At 



,m + l , o+j,m + l 



+ i + /?/ ^ + J1p™J~2 ) (10) 



where 



Ax 4 



4 2 ^-Vj)±-tvj 



(11) 
(12) 



Ax 4 Ax 2 

a = = -Kx~ 4 ~ -h {E ~ ^ + {E ~ 1/3)2 ~ v 'i (13) 

and Vj , Vj and Vj' are the values of the potential and its 
first two space-derivatives at Xj. Due to the boundary 
conditions we can rewrite (10) in a compact form suitable 
for numerical solution 



(14) 



where the matrix 1Z is pentadiagonal with nonvanish- 
ing elements TZa = 1 + At a 8 -, 1Zu±\ = At fif and 
lZa±2 = At 7. Starting with a known ip® the system 
(14) is efficiently solved with standard decomposition and 
back-substitution method [5] in a number of operations 
per time step proportional to J . 

Following the von Neumann stability analysis [5], the 
eigenmodes ^{k) = ^ m e lk J Ax substituted back into (14) 
give 



m + l 



1 



■i6)A*| 



(15) 



where 



a = { y 3 -Ey + —{V 3 -E) s ^ 



kAa 



16 . 2 {kAa 



■ sin 2 (kAx) - V" (16) 



Ax 4 V 2 J Ax 4 y ' 3 
2 

b = V- sin (kAx) . 

Ax 3 y ' 



(17) 



Stability is obtained when the growing ratio (15) is 
smaller than unity, i.e. when 



At > 



2a 



a 2 + b 2 



(18) 



This means that we can choose At very large and obtain 
convergence to the relaxed quantities (5-6) in a single 
iteration of Eq. (14). Differently stated the rate of con- 
vergence for the recursive equation ip m+1 = TZ -1 ^" 1 is 
infinite. This can be seen directly by the formula for the 
rate of convergence [6], — In S^ -1 ), where S(. . .) means 
spectral radius. In the limit At — > 00 the eigenvalues of 
TZT 1 and its spectral radius vanish and the rate of con- 
vergence diverges. Finite computer accuracy imposes a 
limitation on the value of At. The relaxed wavefunction 
is obtained after normalization of the vanishing ip(x,t; E) 
and At cannot be so large that ipj yields underflow. The 
limitation, however, is not crucial and full relaxation can 
be usually obtained with very few iterations (we never 
use more than ten iterations) even without making an 
optimal choice for At (maximum value allowed) . 

Efficiency and precision obtainable from the selective 
relaxation method [7] have been checked in various cases. 
These include exactly solvable problems as well as prob- 
lems where results obtained with different numerical pro- 
cedures are available for comparison [8,9]. Here we report 
only on a comparison with the exact results of the Morse 
potential and on the possibility to evaluate with great 
accuracy the fundamental energy splitting of double-well 
potentials. 

The eigenvalue problem with the Morse potential 



V(x) 



2e" 



(19) 



has well known analytical solutions [10]. In Fig.s 1 and 
2 we compare the exact eigenvalues and eigenfunctions 
with the corresponding energies and wavefunctions ob- 
tained with the selective relaxation method for different 
values of the lattice step Ax. We chose /j, = 0.2 which 
corresponds to have five bound states n = 0, . . . , 4. The 
convergence to the nth bound state is absolutely insen- 
sitive to the choice of the initial wave function as well as 
of the selecting energy E in the interval T n . 

According to the discretization scheme used in (10) 
the algorithm is first-order accurate in the lattice step 
Ax. More explicitly, for the eigenvalues we observe a 
systematic error 



rel 



E n = e n Ax 



(20) 
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with e n increasing with n. Since the computation time 
is proportional to the number of steps in the space lat- 
tice, we conclude that the error in the evaluation of the 
eigenvalues decreases quadratically with the computation 
time. In order to fix the ideas, the computation time for 
J = 10 4 , which is a typical figure in evaluating E4 with a 
0.01 % error, is about 2 s per time iteration in a 25 MHz 
486 PC. 

The increase of e n with n may become a problem when 
evaluating high energy eigenstates. Indeed, the dimen- 
sion J of the space lattice necessary for controlling the 
error (20) through very small Ax may exceed the com- 
puter capacity. The problem is overcome by resorting 
to a higher-order approximation in the discretization of 
the operator (H — E) 2 . If we substitute the right hand 
side of Eq. (10) with a fcth-order accurate expression 
the error (20) gets proportional to Ax k+1 and we have 
higher precision for a given lattice dimension J . In this 
case, however, the matrix 1Z has 2(k + 1) + 1 nonvan- 
ishing diagonals (for x £ R) and the computation time 
for the same J increases [5]. A quantitative comparison 
between the minimal-accuracy algorithm presented here 
and higher-accuracy ones, also in the cases x £ M 2 and 
x £ M 3 is, deferred elsewhere. 

Even in its minimal-accuracy version the selective re- 
laxation method allows us to make a relevant advance 
in the evaluation of the fundamental energy splitting 
T = Ei — Eo of a double- well potential. It is well known 
that this problem gets rapidly unapproachable with stan- 
dard numerical methods when the tunneling rate between 
the two wells decreases [11]. A different situation arises 
for the fundamental energy splitting T re ; obtained as dif- 
ference of the lowest two relaxed eigenvalues 

T rel =T+{e 1 -e )Ax 2 . (21) 

T re i shows a systematic error due to the finite lattice, 
(ei — eo) Ax 2 , which is the difference of two close numbers 
and vanishes for T vanishing. In addition, the condition 
Eo — Ei causes no trouble in selecting the two eigen- 
values since the corresponding eigenfunctions ipo and ipi 
have different parity. Initial wavefunctions ip(x,0;E) 
with the same selecting energy, e.g. chosen as the WKB 
approximation to the ground state, but different parity 
automatically relax toward the eigenfunctions with the 
corresponding parity. 

An example of the discussed behavior is shown in Fig. 
3. We have considered the bistable potential 

V(x) = -Xx 2 + x 4 (22) 

for A = 15. The values of T re ; obtained for different 
values of the space step Ax follows accurately law (21). 
A linear fit gives T = 1.9496 x 10" 10 and e 1 - e = 
2.6404 x 10 -8 . This last figure should be compared with 
the single level accuracy to — £1 — 7.22 which is eight 
orders of magnitude greater. 



Comparison of the splitting values obtained through 
selective relaxation with those obtained through other 
techniques establishes a clear superiority of our method. 
Recently a substantial advance in evaluating the funda- 
mental energy splitting of double- well potentials was real- 
ized with the technique of supersymmetric quantum me- 
chanics [11]. Using as input the ground-state wavefunc- 
tion ipo(%), obtained, for instance, with standard Runge- 
Kutta integration, the supersymmetric approach allows 
one to evaluate the splitting through a logarithmic per- 
turbation series which converges rapidly in the limit of 
small tunneling rate. In Table I we compare the funda- 
mental energy splitting of the double-well potential (22) 
evaluated with the Runge-Kutta method, the supersym- 
metric series at third order and the selective relaxation 
for different values of A. For A small the potential is 
weakly bistable and the Runge-Kutta method is reliable. 
The supersymmetric result shows weak convergence in 
this limit. For A large the Runge-Kutta method is unre- 
liable for evaluating the splitting. However, this method 
is still good for evaluating the ground state used in the su- 
persymmetric series which shows good convergence. The 
selective relaxation method gives the right result in the 
full range of A values with a minimal amount of compu- 
tation time. Notice that, in addition to the cases consid- 
ered in Ref. [11], we are able to evaluate the splitting for 
A = 15, i.e. in a region where T is close to the maximal 
accuracy available in our computer (double precision) . 

Some final comments are in order. We have considered 
only Hamiltonians with nondegenerate spectra. The se- 
lective relaxation method, however, applies also in pres- 
ence of degeneracy which may occur in two- or three- 
dimensional cases. Initial wavefunctions which are mu- 
tually orthogonal converge to the corresponding degen- 
erate eigenfunctions. In fact, this property was used for 
selecting the nearly-degenerate fundamental levels of the 
double- well potential (22). 

Modified boundary conditions, e.g. asymptotic known 
values at the points x m i„ and x max , allow one to evaluate 
eigenfunctions in the continuous spectrum. In particular, 
extended resonance states can be found through the so 
called quantum transmitting boundary method [12]. Pe- 
riodic boundary conditions allow one to calculate eigen- 
values and eigenfunctions of periodic potentials. 

The main characteristics of the proposed relaxation 
method, namely selectivity and absolute stability, can be 
extended to different classes of eigenvalue problems such 
as those involving Fokker-Planck and Dirac operators. 
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TABLE I. Fundamental energy splitting of the double- well 
potential V(x) = — \x 2 + x 4 obtained in Ref. [11] with 
Runge-Kutta integration, Trk, and supersymmetric pertur- 
bation series at third order, Tss3 , compared with the selective 
relaxation result, T re i, at lattice step Ax = 10 -3 . Notice that 
numerical Runge-Kutta calculations are unreliable for A > 10 
while supersymmetric result gets inaccurate for A small. 
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FIG. 1. Comparison of the exact eigenvalues E„ of the 
Morse potential with the relaxed eigenvalues E rF j for dif- 
ferent values of the lattice step Ax. The potential is 
V(x) = e _2MX - 2e _MX with /j = 0.2 and has five bound states 
n = 0, . . . ,4. 
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FIG. 2. Comparison of the exact eigenfunctions tp n of the 
potential of Fig. 1 with the relaxed eigenfunctions tp rF j for 
different values of the lattice step Ax. 
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FIG. 3. Fundamental energy splitting T rF j obtained as 
difference of the lowest two relaxed eigenvalues of the 
potential V(x) = —15a; 2 + x 4 for different values of 
the lattice step Ax. A linear fit (solid line) gives 
T re i = 1.9496 x 10" 10 + 2.6404 x 10" 8 Ax 2 . 
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